library(ggplot2)
library(ggpubr)
theme_customize=theme_bw()+
  theme(panel.border=element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.line = element_line(colour = "black"),
        legend.title = element_blank())+
  theme(axis.text.x=element_text(size=14,angle=0,color="Black"),
        axis.text.y=element_text(size=14,angle=0,color="Black"))+
  theme(axis.title.x = element_text(size = 15, color = "Black"),
        axis.title.y = element_text(size = 15,  color = "Black"))

multiplot <- function(..., plotlist=NULL, file, cols=1, layout=NULL) {
  library(grid)
  
  # Make a list from the ... arguments and plotlist
  plots <- c(list(...), plotlist)
  
  numPlots = length(plots)
  
  # If layout is NULL, then use 'cols' to determine layout
  if (is.null(layout)) {
    # Make the panel
    # ncol: Number of columns of plots
    # nrow: Number of rows needed, calculated from # of cols
    layout <- matrix(seq(1, cols * ceiling(numPlots/cols)),
                     ncol = cols, nrow = ceiling(numPlots/cols))
  }
  
  if (numPlots==1) {
    print(plots[[1]])
    
  } else {
    # Set up the page
    grid.newpage()
    pushViewport(viewport(layout = grid.layout(nrow(layout), ncol(layout))))
    
    # Make each plot, in the correct location
    for (i in 1:numPlots) {
      # Get the i,j matrix positions of the regions that contain this subplot
      matchidx <- as.data.frame(which(layout == i, arr.ind = TRUE))
      
      print(plots[[i]], vp = viewport(layout.pos.row = matchidx$row,
                                      layout.pos.col = matchidx$col))
    }
  }
}

###bar graph of trophic groups
setwd("D:/入侵/数据/分析2021/4 trophic groups")
##variance analysis 
library(vegan)
library("reshape2")
trophic_den=read.csv("trophic_den.csv",header=TRUE)
dat1<-melt(trophic_den,
           id.vars = c('Community','Site'),
           variable.name='ind',
           value.name='val')

library(dplyr)
list_name1 <-  unique(dat1$ind)
##normality    ########None of them are normal     ####>0.05 normality
sink("normal-test result.txt")  #Save the output
for(i in list_name1){
  print(paste(i))
  tmp<-shapiro.test((dat1 %>% filter(ind == i))$val)
  print(tmp)
}
sink()

library(agricolae)
sink("Kruskal result.txt")  #Save the output
for(i in list_name1){
  print(paste(i))
  mod<-dat1 %>% filter(ind == i) %>% with(.,kruskal(val,Community,group=TRUE, main=paste(i,sep = " ")))
  print(mod)
}
sink()

bar_trop=read.csv("bar_trop.csv",header=TRUE)
list_trop <-  unique(bar_trop$ind)
for(i in list_trop){
  tmp<-bar_trop %>% 
    filter(ind == i) %>% 
    dcast(., Community ~ type, value.var = "val") 
  tmp$Community = factor(tmp$Community, levels=c('OP','IP','IPS','IS','RP'))
  tmp$mean=as.numeric(tmp$mean)
  tmp$SE=as.numeric(tmp$SE)
  assign(paste("p_num_", i, sep = ""), ggplot(tmp,aes(x=Community,y=mean,fill=Community))+
           theme_customize+
           theme(legend.position="none")+
           geom_bar(position="dodge",stat="identity",width=0.5)+
           geom_text(aes(y=mean+3*SE,label=text),position=position_dodge(width = 0.9))+ #柱状图上添加文字
           xlab("Plant Community")+
           ylab(paste(i))+
           geom_errorbar(aes(ymax = mean+1.96*SE, ymin = mean-1.96*SE), 
                         position = position_dodge(0.9), width = 0.15))
}
for(i in list_trop){
  assign(paste("p_num_",i, sep = ""),
         get(paste("p_num_",i, sep = ""))+coord_cartesian(ylim = c(0,40)))
}
multiplot(p_num_detritivores, p_num_predators,  
          p_num_chewers, p_num_parasites,
          p_num_sap_feeders,  p_num_web_builders,
          p_num_stem_feeders, p_num_hunters, cols=4)  ##16*8